#############################################
##### Replication Code
##### Mobilizing for Prisoner Release During a Pandemic via Perspective-Taking Exercises: Appendix 3
##### Mackenzie Israel-Trummel
#############################################

##First run the replication file for the Manuscript 


###############
##Survey Completion Time in Appendix 1
###############

summary(covid.white$Duration..in.seconds.)/60


###############
##Sample Descriptives
###############

##Clean demographic variables and all other variables used in this file

table(covid.white$agecat)
covid.white$age[covid.white$agecat == "18_24"] <- 0
covid.white$age[covid.white$agecat == "25_34"] <- 1
covid.white$age[covid.white$agecat == "35_44"] <- 2
covid.white$age[covid.white$agecat == "45_54"] <- 3
covid.white$age[covid.white$agecat == "55_64"] <- 4
covid.white$age[covid.white$agecat == "65+"] <- 5
table(covid.white$agecat, covid.white$age)
covid.white$age <- as.numeric(covid.white$age)

table(covid.white$Q91)
covid.white$woman[covid.white$Q91 == "Man"] <- 0
covid.white$woman[covid.white$Q91 == "Woman"] <- 1
table(covid.white$Q91, covid.white$woman)

table(covid.white$Q9)
covid.white$edu[covid.white$Q9 == "Grades 1-8"] <- 0
covid.white$edu[covid.white$Q9 == "Some high school"] <- 0
covid.white$edu[covid.white$Q9 == "High school graduate or GED"] <- 1
covid.white$edu[covid.white$Q9 == "Some college, or 2-year degree"] <- 2
covid.white$edu[covid.white$Q9 == "4-year college graduate"] <- 3
covid.white$edu[covid.white$Q9 == "Post-graduate education"] <- 4
table(covid.white$Q9, covid.white$edu)


table(covid.white$Q64_1)
covid.white$resent_deserve[covid.white$Q64_1 == "Strongly agree"] <- 0
covid.white$resent_deserve[covid.white$Q64_1 == "Somewhat agree"] <- 1
covid.white$resent_deserve[covid.white$Q64_1 == "Neither agree nor disagree"] <- 2
covid.white$resent_deserve[covid.white$Q64_1 == "Somewhat disagree"] <- 3
covid.white$resent_deserve[covid.white$Q64_1 == "Strongly disagree"] <- 4
table(covid.white$Q64_1, covid.white$resent_deserve)

table(covid.white$Q64_2)
covid.white$resent_nofavors[covid.white$Q64_2 == "Strongly agree"] <- 4
covid.white$resent_nofavors[covid.white$Q64_2 == "Somewhat agree"] <- 3
covid.white$resent_nofavors[covid.white$Q64_2 == "Neither agree nor disagree"] <- 2
covid.white$resent_nofavors[covid.white$Q64_2 == "Somewhat disagree"] <- 1
covid.white$resent_nofavors[covid.white$Q64_2 == "Strongly disagree"] <- 0
table(covid.white$Q64_2, covid.white$resent_nofavors)

table(covid.white$Q64_3)
covid.white$resent_tryharder[covid.white$Q64_3 == "Strongly agree"] <- 4
covid.white$resent_tryharder[covid.white$Q64_3 == "Somewhat agree"] <- 3
covid.white$resent_tryharder[covid.white$Q64_3 == "Neither agree nor disagree"] <- 2
covid.white$resent_tryharder[covid.white$Q64_3 == "Somewhat disagree"] <- 1
covid.white$resent_tryharder[covid.white$Q64_3 == "Strongly disagree"] <- 0
table(covid.white$Q64_3, covid.white$resent_tryharder)


table(covid.white$Q64_4)
covid.white$resent_legacy[covid.white$Q64_4 == "Strongly agree"] <- 0
covid.white$resent_legacy[covid.white$Q64_4 == "Somewhat agree"] <- 1
covid.white$resent_legacy[covid.white$Q64_4 == "Neither agree nor disagree"] <- 2
covid.white$resent_legacy[covid.white$Q64_4 == "Somewhat disagree"] <- 3
covid.white$resent_legacy[covid.white$Q64_4 == "Strongly disagree"] <- 4
table(covid.white$Q64_4, covid.white$resent_legacy)

covid.white$resentment <- covid.white$resent_deserve + covid.white$resent_nofavors + covid.white$resent_tryharder + covid.white$resent_legacy
table(covid.white$resentment)
summary(covid.white$resentment)
quantile(covid.white$resentment, probs = seq (from = .33, to = .66, by =.33), na.rm = TRUE)
covid.white$resentment3[covid.white$resentment < 7] <- 0
covid.white$resentment3[covid.white$resentment > 6 &  covid.white$resentment < 10] <- 1
covid.white$resentment3[covid.white$resentment > 9] <- 2
table(covid.white$resentment, covid.white$resentment3)




table(covid.white$Q53_3)
covid.white$poordeserve[covid.white$Q53_3 == "More money than they deserve"] <- 0
covid.white$poordeserve[covid.white$Q53_3 == "About the right amount of money"] <- 1
covid.white$poordeserve[covid.white$Q53_3 == "Less money than they deserve"] <- 2
table(covid.white$Q53_3, covid.white$poordeserve)

covid.white$poordeservemore[covid.white$Q53_3 == "More money than they deserve"] <- 0
covid.white$poordeservemore[covid.white$Q53_3 == "About the right amount of money"] <- 0
covid.white$poordeservemore[covid.white$Q53_3 == "Less money than they deserve"] <- 1
table(covid.white$Q53_3, covid.white$poordeservemore)


table(covid.white$Q65)
covid.white$whitefate[covid.white$Q65 == "No"] <- 0
covid.white$whitefate[covid.white$Q65 == "Yes, but not very much"] <- 1
covid.white$whitefate[covid.white$Q65 == "Yes, some"] <- 2
covid.white$whitefate[covid.white$Q65 == "Yes, a lot"] <- 3
table(covid.white$Q65, covid.white$whitefate)

table(covid.white$Q68)
covid.white$covidworry[covid.white$Q68 == "Not at all worried"] <- 0
covid.white$covidworry[covid.white$Q68 == "A little worried"] <- 1
covid.white$covidworry[covid.white$Q68 == "Pretty worried"] <- 2
covid.white$covidworry[covid.white$Q68 == "Very worried"] <- 3
table(covid.white$Q68, covid.white$covidworry)



###############
## Table A1: Sample Descriptives
###############

stargazer(covid.white[c("woman", "age", "edu", "pid7")], omit.summary.stat = c("sd")) 
##add proportions for region variable manually
table(covid.white$region)
table(covid.white$region)/1534



###############
## Table A2: Experimental Treatments Balance Table
###############

##Produce table: means across conditions

control <- subset(covid.white, condition == "control")
notreat <- subset(covid.white, condition == "notreat")
t1 <- subset(covid.white, condition == "t1")
t2 <- subset(covid.white, condition == "t2")


v1 <- notreat[c("age", "woman", "edu", "resentment", "poordeserve", "pid7", "covidworry", "whitefate")]
v2 <- control[c("age", "woman", "edu", "resentment", "poordeserve", "pid7", "covidworry", "whitefate")]
v3 <- t1[c("age", "woman", "edu", "resentment", "poordeserve", "pid7", "covidworry", "whitefate")]
v4 <- t2[c("age", "woman", "edu", "resentment", "poordeserve", "pid7", "covidworry", "whitefate")]

stargazer(v1, omit.summary.stat = c("n" , "sd", "min", "max")) 
stargazer(v2, omit.summary.stat = c("n" , "sd", "min", "max")) 
stargazer(v3, omit.summary.stat = c("n" , "sd", "min", "max")) 
stargazer(v4, omit.summary.stat = c("n" , "sd", "min", "max")) 


chisq.test(covid.white$age, covid.white$condition)
chisq.test(covid.white$woman, covid.white$condition)
chisq.test(covid.white$edu, covid.white$condition)
chisq.test(covid.white$resentment, covid.white$condition)
chisq.test(covid.white$poordeserve, covid.white$condition)
chisq.test(covid.white$pid7, covid.white$condition)
chisq.test(covid.white$covidworry, covid.white$condition)
chisq.test(covid.white$whitefate, covid.white$condition)
table(covid.white$condition)
chisq.test(table(covid.white$condition))






###############
##Emotional Responses to the Perspective-Taking treatments as compared to the Control condition
###############

table(covid$Q51) ##control feelings
table(covid$Q52) ##egocentric feelings
table(covid$Q53) ##surrogate feelings

load(file = "emotion.rdata")

emotions <- merge(covid, emotion, by = "id")

aggregate(concern ~ condition, data = emotions, FUN = mean)
emotion.mod <- lm(concern ~ condition, data = emotions)
summary(emotion.mod)

emotions.white <- subset(emotions, white == 1)
aggregate(concern ~ condition, data = emotions.white, FUN = mean)
emotion.mod.w <- lm(concern ~ condition, data = emotions.white)
summary(emotion.mod.w)


###########
#####Table A3: Emotional Responses by Treatment Condition
###########

texreg(list(emotion.mod, emotion.mod.w), 
       single.row = TRUE, 
       custom.coef.names = c("(Intercept)", 
                             "Perspective: Egocentric", 
                             "Perspective: Surrogate"
       ), 
       stars = 0.05, dcolumn=TRUE, caption = "Emotional Responses By Treatment Condition", label = "tab:emotions",
       custom.model.names = c("All Respondents", "White Respondents"), 
       booktabs = TRUE)


###########
##Table A4: Treatment Means
##########

aggregate(covid.white$release, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE) 
aggregate(covid.white$supportrelease_Rec3, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)

##corrections for multiple comparisons
pairwise.t.test(covid.white$release, covid.white$condition, p.adjust.method="holm")
pairwise.t.test(covid.white$supportrelease_Rec3, covid.white$condition, p.adjust.method="holm")


########
## Modeling Treatment Effects on Individual Items in the Policy Opinion Scale
#######

##individual items

mod.all <- lm(release_all ~ condition, data = covid.white)
summary(mod.all)

mod.old <- lm(release_old ~ condition, data = covid.white)
summary(mod.old)

mod.immuno <- lm(release_immuno ~ condition, data = covid.white)
summary(mod.immuno)

mod.nonviolent <- lm(release_nonviolent ~ condition, data = covid.white)
summary(mod.nonviolent)

mod.bail <- lm(release_bail ~ condition, data = covid.white)
summary(mod.bail)


###########
##Table A5: Treatment Effects on Individual Policy Opinion Items
##########


texreg(list(mod.bail, mod.old, mod.immuno, mod.nonviolent, mod.all), 
       single.row = TRUE, 
       custom.coef.names = c("(Intercept)", 
                             "Informational Control",
                             "Perspective: Egocentric", 
                             "Perspective: Surrogate"
       ), 
       stars = 0.05, dcolumn = TRUE,
       custom.model.names = c("Bail", "Elderly", "Immuno.", "Nonviolent", "All prisoners"), caption = c("Treatment Effects on Individual Release Attitudes Items"), label = c("tab:items"), 
       booktabs = TRUE)




##############
##Binned Release Attitudes
##############
table(covid.white$release)
quantile(covid.white$release, probs = seq (from = 0, to = 1, by =.20), na.rm = TRUE)
covid.white$release_binned[covid.white$release == 0] <- 0
covid.white$release_binned[covid.white$release > 0 & covid.white$release < 6.4] <- 1
covid.white$release_binned[covid.white$release > 6.4 & covid.white$release < 10] <- 2
covid.white$release_binned[covid.white$release > 9 & covid.white$release < 14] <- 3
covid.white$release_binned[covid.white$release > 13 ] <- 4
table(covid.white$release, covid.white$release_binned)
table(covid.white$release_binned)

aggregate(covid.white$release_binned, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE) 
aggregate(covid.white$release_binned, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[2] - aggregate(covid.white$release_binned, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[1, 2] ##diff from no treat

aggregate(covid.white$release_binned, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[2] - aggregate(covid.white$release_binned, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[2, 2] ##diff from control


## subsets by treatment
control <- subset(covid.white, condition == "control")
notreat <- subset(covid.white, condition == "notreat")
t1 <- subset(covid.white, condition == "t1")
t2 <- subset(covid.white, condition == "t2")


#t-tests between treatments, notreat as baseline
t.a <- t.test(control$release_binned, notreat$release_binned)
t.b <- t.test(t1$release_binned, notreat$release_binned)
t.c <- t.test(t2$release_binned, notreat$release_binned) 

t.test(t1$release_binned, control$release_binned)
t.test(t2$release_binned, control$release_binned)
t.test(t2$release_binned, t1$release_binned)



pairwise.t.test(covid.white$release_binned, covid.white$condition, p.adjust.method="holm")


treateffect.release <- c(0, t.a$estimate[1]- t.a$estimate[2], t.b$estimate[1]- t.b$estimate[2], t.c$estimate[1]- t.c$estimate[2])

lower.release <- c(0, t.a$conf.int[1], t.b$conf.int[1], t.c$conf.int[1])

upper.release <- c(0, t.a$conf.int[2], t.b$conf.int[2], t.c$conf.int[2])


treatment <- c("No\n treatment", "Informational\n Control",  "Perspective:\n egocentric", "Perspective:\n surrogate")

data.release = as.data.frame(treateffect.release)
data.release$lower = lower.release
data.release$upper = upper.release
data.release$treatment = treatment
data.release$treatment = factor(data.release$treatment, levels = c("Perspective:\n surrogate", "Perspective:\n egocentric",  "Informational\n Control", "No\n treatment"))
data.release$labels <- round(data.release$treateffect.release, 3)
data.release$labels[data.release$labels == 0] <- NA

p_append1 <- ggplot(data = data.release, aes(x = treateffect.release, y = treatment)) +
  geom_text(aes(label = labels), nudge_y = 0.2, nudge_x = -0.002, size=6) +
  theme_bw() +
  geom_point(size = 4) +
  geom_errorbarh(aes(y = treatment, xmin = lower, xmax = upper, height=.1),size=.5) +
  xlab("Effect on Support for Releasing Prisoners")  + ylab("") +
  theme(axis.text = element_text(size = 18), #Increase size of axis text
        axis.title = element_text(size = 18))	+  geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") +
  #ggtitle(title)+
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(0.1, .5, 0.1, 0), "cm")) +
  annotate("text", x = 0.1, y = 4.15, parse = TRUE, label = "bar(x)==1.950", col = "black", size = 6) 


pdf("figA1.pdf", height= 5, width = 7)
p_append1
dev.off()





##############
##Alternative models for 3 point behavioral outcome
##############


beh_olr <- polr(as.factor(supportrelease_Rec3) ~ condition, data = covid.white, Hess= T)
summary(beh_olr)


beh_lm <- lm(supportrelease_Rec3 ~ condition, data = covid.white)
summary(beh_lm)



###########
##Table A6: Alternative Modeling of Semi-Behavioral Effects
##########

texreg(list(beh_lm, beh_olr),
       single.row = TRUE, 
       custom.coef.names = c("(Intercept)", 
                             "Informational Control",
                             "Perspective: Egocentric", 
                             "Perspective: Surrogate",
                             "Cut 1", 
                             "Cut 2"
       ), 
       stars = 0.05, dcolumn = TRUE, 
       custom.model.names = c("OLS", "Ordered Logit"), caption = c("Alternative Modeling of Semi-Behavioral Effects"), label = c("tab:robust_semibeh"), 
       booktabs = TRUE)




################
## Alternative Codings of Semibehavioral
###############

##Modal Coding

table(covid.white$Disagree, covid.white$Modal) ##Modal score among the ones with any disagreement

##Add the agreed upon scores to the Modal variable
covid.white$Modal[covid.white$Disagree == "Agree" & covid.white$Mackenzie == 1] <- 1
covid.white$Modal[covid.white$Disagree == "Agree" & covid.white$Mackenzie == 2] <- 2
covid.white$Modal[covid.white$Disagree == "Agree" & covid.white$Mackenzie == 3] <- 3
covid.white$Modal[covid.white$Disagree == "Agree" & covid.white$Mackenzie == 4] <- 4
covid.white$Modal[covid.white$Disagree == "Agree" & covid.white$Mackenzie == 5] <- 5
covid.white$Modal[covid.white$Disagree == "Agree" & covid.white$Mackenzie == 99] <- 99
table(covid.white$Modal)
table(covid.white$Disagree, covid.white$Modal)

covid.white$supportrelease_Modal3 <- 0
covid.white$supportrelease_Modal3[covid.white$Modal == 1] <- 1
covid.white$supportrelease_Modal3[covid.white$Modal == 2] <- -1
table(covid.white$Modal, covid.white$supportrelease_Modal3)

##Single coder
table(covid.white$Mackenzie)

covid.white$singlecoder3 <- 0
covid.white$singlecoder3[covid.white$Mackenzie == 1] <- 1
covid.white$singlecoder3[covid.white$Mackenzie == 2] <- -1
table(covid.white$Mackenzie, covid.white$singlecoder3)


############
##Table A7: Alternative Codings of Semibheavioral DV: Means by Treatment Group
############

aggregate(covid.white$supportrelease_Modal3, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)
pairwise.t.test(covid.white$supportrelease_Modal3, covid.white$condition, p.adjust.method="holm")

aggregate(covid.white$singlecoder3, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE) 
pairwise.t.test(covid.white$singlecoder3, covid.white$condition, p.adjust.method="bonferroni")




###########
## Means by Moderator Groups 
###########


###########
## Table A8: Conditional Treatment Means: Carceral Contact
###########

##Carceral contact
aggregate(contact$release, by = list(Category = contact$condition), FUN = mean, na.rm = TRUE)
aggregate(contact$supportrelease_Rec3, by = list(Category = contact$condition), FUN = mean, na.rm = TRUE)


##No carceral contact
aggregate(nocontact$release, by = list(Category = nocontact$condition), FUN = mean, na.rm = TRUE) 
aggregate(nocontact$supportrelease_Rec3, by = list(Category = nocontact$condition), FUN = mean, na.rm = TRUE) 

##t-tests
contact.nt <- subset(contact, condition == "notreat")
contact.c <- subset(contact, condition == "control")
contact.t1 <- subset(contact, condition == "t1")
contact.t2 <- subset(contact, condition == "t2")

nocontact.nt <- subset(nocontact, condition == "notreat")
nocontact.c <- subset(nocontact, condition == "control")
nocontact.t1 <- subset(nocontact, condition == "t1")
nocontact.t2 <- subset(nocontact, condition == "t2")

### Contact: opinion
##compared to no treatment group 
t.test(contact.nt$release, contact.c$release)
t.test(contact.nt$release, contact.t1$release)
t.test(contact.nt$release, contact.t2$release)

##compared to control group
t.test(contact.c$release, contact.t1$release)
t.test(contact.c$release, contact.t2$release)


### Contact: behavioral
##compared to no treatment group 
t.test(contact.nt$supportrelease_Rec3, contact.c$supportrelease_Rec3)
t.test(contact.nt$supportrelease_Rec3, contact.t1$supportrelease_Rec3)
t.test(contact.nt$supportrelease_Rec3, contact.t2$supportrelease_Rec3)

##compared to control group
t.test(contact.c$supportrelease_Rec3, contact.t1$supportrelease_Rec3)
t.test(contact.c$supportrelease_Rec3, contact.t2$supportrelease_Rec3)


### No Contact: opinion
##compared to no treatment group 
t.test(nocontact.nt$release, nocontact.c$release)
t.test(nocontact.nt$release, nocontact.t1$release)
t.test(nocontact.nt$release, nocontact.t2$release)

##compared to control group
t.test(nocontact.c$release, nocontact.t1$release)
t.test(nocontact.c$release, nocontact.t2$release)


### No Contact: behavioral
##compared to no treatment group 
t.test(nocontact.nt$supportrelease_Rec3, nocontact.c$supportrelease_Rec3)
t.test(nocontact.nt$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)
t.test(nocontact.nt$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)

##compared to control group
t.test(nocontact.c$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)
t.test(nocontact.c$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)




##comparing treated no contact groups to contact groups for semibehavioral measure



t.test(contact.nt$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)
t.test(contact.nt$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)

t.test(contact.c$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)
t.test(contact.c$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)

t.test(contact.t1$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)
t.test(contact.t1$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)
t.test(contact.t2$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)
t.test(contact.t2$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)



###########
## Table A9: Conditional Treatment Means: Partisanship
###########


##Democrats
aggregate(dems$release, by = list(Category = dems$condition), FUN = mean, na.rm = TRUE)
aggregate(dems$supportrelease_Rec3, by = list(Category = dems$condition), FUN = mean, na.rm = TRUE)


##Republicans
aggregate(reps$release, by = list(Category = reps$condition), FUN = mean, na.rm = TRUE) 
aggregate(reps$supportrelease_Rec3, by = list(Category = reps$condition), FUN = mean, na.rm = TRUE) 


dems.nt <- subset(dems, condition == "notreat")
dems.c <- subset(dems, condition == "control")
dems.t1 <- subset(dems, condition == "t1")
dems.t2 <- subset(dems, condition == "t2")

t.test(dems.nt$release, dems.c$release)
t.test(dems.nt$release, dems.t1$release)
t.test(dems.nt$release, dems.t2$release)

t.test(dems.c$release, dems.t1$release)
t.test(dems.c$release, dems.t2$release)

t.test(dems.nt$supportrelease_Rec3, dems.c$supportrelease_Rec3)
t.test(dems.nt$supportrelease_Rec3, dems.t1$supportrelease_Rec3)
t.test(dems.nt$supportrelease_Rec3, dems.t2$supportrelease_Rec3)

t.test(dems.c$supportrelease_Rec3, dems.t1$supportrelease_Rec3)
t.test(dems.c$supportrelease_Rec3, dems.t2$supportrelease_Rec3)



reps.nt <- subset(reps, condition == "notreat")
reps.c <- subset(reps, condition == "control")
reps.t1 <- subset(reps, condition == "t1")
reps.t2 <- subset(reps, condition == "t2")

t.test(reps.nt$release, reps.c$release)
t.test(reps.nt$release, reps.t1$release)
t.test(reps.nt$release, reps.t2$release)

t.test(reps.c$release, reps.t1$release)
t.test(reps.c$release, reps.t2$release)

t.test(reps.nt$supportrelease_Rec3, reps.c$supportrelease_Rec3)
t.test(reps.nt$supportrelease_Rec3, reps.t1$supportrelease_Rec3)
t.test(reps.nt$supportrelease_Rec3, reps.t2$supportrelease_Rec3)

t.test(reps.c$supportrelease_Rec3, reps.t1$supportrelease_Rec3)
t.test(reps.c$supportrelease_Rec3, reps.t2$supportrelease_Rec3)



##########
## Table A10: Regression Models for CATEs 
##########

summary(contact_att_lm)
summary(contact_beh_lm)
summary(pid_att_lm)
summary(pid_beh_lm)



texreg(list(contact_att_lm, contact_beh_lm, pid_att_lm, pid_beh_lm),
       single.row = FALSE, 
       custom.coef.names = c("(Intercept)", 
                             "Informational Control",
                             "Perspective: Egocentric", 
                             "Perspective: Surrogate",
                             "Contact", 
                             "Inf. Control*Contact",
                             "Egocentric*Contact",
                             "Surrogate*Contact",
                             "Democrat", 
                             "Inf. Control*Democrat",
                             "Egocentric*Democrat",
                             "Surrogate*Democrat"
       ), 
       stars = 0.05, 
       custom.model.names = c("Policy", "Semibehavioral", "Policy", "Semibehavioral"), caption = c("Regression Models for CATEs"), label = c("tab:cates"), 
       booktabs = TRUE, dcolumn = TRUE)



################
## Figure A2: Conditional Treatment Effects by Racial Resentment
###############



resent_att_lm <- lm(release ~ condition*resentment3, data = covid.white)
summary(resent_att_lm)

ef_resent_att <- effect("condition : resentment3", resent_att_lm)

ef_resent_att

d_resent_att <- as.data.frame(ef_resent_att) #save as data frame
d_resent_att

#Save only estimates at 0, 1, 2
d_resent_att <- subset(d_resent_att, resentment3 == 0 | resentment3 == 1 | resentment3 == 2)
d_resent_att



#Add in factor names for plotting
d_resent_att$names <- recode(d_resent_att$resentment3, '0' = 'Low Resentment' , '1' = 'Medium Resentment', '2' = 'High Resentment')
d_resent_att
d_resent_att$names<- factor(d_resent_att$names, levels = c("Low Resentment", "Medium Resentment", "High Resentment"))

d_resent_att$cond_lab <- recode(d_resent_att$condition, 'notreat' = 'No\n Treatment' , 'control' = 'Informational\n Control', 't1' = 'Perspective:\n Egocentric', 't2' = 'Perspective:\n Surrogate')
d_resent_att


#Plot and Format
resentplot <- ggplot(d_resent_att, aes(x = cond_lab, y = fit, shape = names, color = names)) +  
  scale_color_manual(values=c("#ffd353", "#ef8737", "#bb292c"), name=" ") + 
  geom_point(size = 4, position = position_dodge(.2)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.15, position = position_dodge(.2)) +
  theme_bw()  +
  ylab("Policy Opinion") + xlab(" ")+
  ylim(4, 13) + 
  scale_shape_manual(values = c(16, 17, 18), name=" ")+
  theme(legend.position = "bottom", legend.text = element_text(size = 16),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 16), 
        axis.text.y = element_text(size = 18), 
        axis.title.y = element_text(size = 18))


pdf("figA2a.pdf", height= 7, width = 7)
resentplot
dev.off()


resent_beh_lm <- lm(supportrelease_Rec3 ~ condition*resentment3, data = covid.white)
summary(resent_beh_lm)

ef_resent_beh <- effect("condition : resentment3", resent_beh_lm)

ef_resent_beh

d_resent_beh <- as.data.frame(ef_resent_beh) #save as data frame
d_resent_beh

#Save only estimates at 1 and 0
d_resent_beh <- subset(d_resent_beh, resentment3 == 0 | resentment3 == 1 | resentment3 == 2)
d_resent_beh



#Add in factor names for plotting
d_resent_beh$names <- recode(d_resent_beh$resentment3, '0' = 'Low Resentment' , '1' = 'Medium Resentment', '2' = 'High Resentment')
d_resent_beh
d_resent_beh$names<- factor(d_resent_beh$names, levels = c("Low Resentment", "Medium Resentment", "High Resentment"))

d_resent_beh$cond_lab <- recode(d_resent_beh$condition, 'notreat' = 'No\n Treatment' , 'control' = 'Informational\n Control', 't1' = 'Perspective:\n Egocentric', 't2' = 'Perspective:\n Surrogate')
d_resent_beh


#Plot and Format
resentplot_beh <- ggplot(d_resent_beh, aes(x = cond_lab, y = fit, shape = names, color = names)) +  
  scale_color_manual(values=c("#ffd353", "#ef8737", "#bb292c"), name=" ") + 
  geom_point(size = 4, position = position_dodge(.2)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.15, position = position_dodge(.2)) + 
  theme_bw()  +
  ylab("Semibehavioral Outcome") + xlab(" ")+
  ylim(-.4, .4) + 
  scale_shape_manual(values = c(16, 17, 18), name=" ")+
  theme(legend.position = "bottom", legend.text = element_text(size = 16),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 16), 
        axis.text.y = element_text(size = 18), 
        axis.title.y = element_text(size = 18))


pdf("FigA2b.pdf", height= 7, width = 7)
resentplot_beh
dev.off()


####################
## Figure A3: Conditional Treatment Effects by Covid Worry
####################



covwor_att_lm <- lm(release ~ condition*covidworry, data = covid.white)
summary(covwor_att_lm)

ef_covwor_att <- effect("condition : covidworry", covwor_att_lm, xlevels=4)
ef_covwor_att

d_covwor_att <- as.data.frame(ef_covwor_att) #save as data frame
d_covwor_att



#Add in factor names for plotting
d_covwor_att$names <- recode(d_covwor_att$covidworry, '0' = 'Not at all' , '1' = 'A little', '2' = 'Pretty', '3' = 'Very')
d_covwor_att
d_covwor_att$names<- factor(d_covwor_att$names, levels = c('Not at all', 'A little', 'Pretty', 'Very'))

d_covwor_att$cond_lab <- recode(d_covwor_att$condition, 'notreat' = 'No\n Treatment' , 'control' = 'Informational\n Control', 't1' = 'Perspective:\n Egocentric', 't2' = 'Perspective:\n Surrogate')
d_covwor_att


#Plot and Format
covworplot <- ggplot(d_covwor_att, aes(x = cond_lab, y = fit, shape = names, color = names)) +  
  scale_color_manual(values=c("#ffd353", "#ef8737", "#bb292c", "#62205f"), name=" ") + 
  geom_point(size = 4, position = position_dodge(.2)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.15, position = position_dodge(.2)) +
  theme_bw()  +
  ylab("Policy Opinion") + xlab(" ")+
  ylim(4, 13) + 
  scale_shape_manual(values = c(16, 17, 18, 15), name=" ")+
  theme(legend.position = "bottom", legend.text = element_text(size = 16),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 16), 
        axis.text.y = element_text(size = 18), 
        axis.title.y = element_text(size = 18))


pdf("FigA3a.pdf", height= 7, width = 7)
covworplot
dev.off()


covwor_beh_lm <- lm(supportrelease_Rec3 ~ condition*covidworry, data = covid.white)
summary(covwor_beh_lm)

ef_covworry_beh <- effect("condition : covidworry", covwor_beh_lm, xlevels=4)

ef_covworry_beh

d_covworry_beh <- as.data.frame(ef_covworry_beh) #save as data frame
d_covworry_beh



#Add in factor names for plotting
d_covworry_beh$names <- recode(d_covworry_beh$covidworry, '0' = 'Not at all' , '1' = 'A little', '2' = 'Pretty', '3' = 'Very')
d_covworry_beh
d_covworry_beh$names<- factor(d_covworry_beh$names, levels = c('Not at all', 'A little', 'Pretty', 'Very'))

d_covworry_beh$cond_lab <- recode(d_covworry_beh$condition, 'notreat' = 'No\n Treatment' , 'control' = 'Informational\n Control', 't1' = 'Perspective:\n Egocentric', 't2' = 'Perspective:\n Surrogate')
d_covworry_beh


#Plot and Format
covworplot_beh <- ggplot(d_covworry_beh, aes(x = cond_lab, y = fit, shape = names, color = names)) +  
  scale_color_manual(values=c("#ffd353", "#ef8737", "#bb292c", "#62205f"), name=" ") + 
  geom_point(size = 4, position = position_dodge(.2)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.15, position = position_dodge(.2)) + 
  theme_bw()  +
  ylab("Semibehavioral Outcome") + xlab(" ")+
  ylim(-.4, .4) + 
  scale_shape_manual(values = c(16, 17, 18, 15), name=" ")+
  theme(legend.position = "bottom", legend.text = element_text(size = 16),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 16), 
        axis.text.y = element_text(size = 18), 
        axis.title.y = element_text(size = 18))


pdf("FigA3b.pdf", height= 7, width = 7)
covworplot_beh
dev.off()


###########
##Table A11: Treatment Effects Among Respondents of Color
###########

covid.nonwhite <- subset(covid, white != 1)
covid.nonwhite$condition <- as.factor(covid.nonwhite$condition)
covid.nonwhite$condition <- relevel(covid.nonwhite$condition, ref = "notreat")

table(covid.nonwhite$white)
table(covid.nonwhite$condition)

aggregate(covid.nonwhite$release, by = list(Category = covid.nonwhite$condition), FUN = mean, na.rm = TRUE) 
pairwise.t.test(covid.nonwhite$release, covid.nonwhite$condition, p.adjust.method="holm")

aggregate(covid.nonwhite$supportrelease_Rec3, by = list(Category = covid.nonwhite$condition), FUN = mean, na.rm = TRUE) 
pairwise.t.test(covid.nonwhite$supportrelease_Rec3, covid.nonwhite$condition, p.adjust.method="holm")


#################
## Table A12: Observational Analysis of Predictors of Semibehavioral Outcome
#################


##rescale all predictors 0 to 1
table(covid.white$age)
covid.white$age <- as.numeric(covid.white$age)
covid.white$age_1 <- covid.white$age/5
table(covid.white$age_1)

table(covid.white$woman)

table(covid.white$edu)
covid.white$edu_1 <- covid.white$edu/4
table(covid.white$edu_1)

table(covid.white$pid7)
covid.white$pid7_1 <- covid.white$pid7/6
table(covid.white$pid7_1)

table(covid.white$covidworry)
covid.white$covidworry_1 <- covid.white$covidworry/3
table(covid.white$covidworry_1)

table(covid.white$poordeserve)
covid.white$poordeserve_1 <- covid.white$poordeserve/2
table(covid.white$poordeserve_1)

table(covid.white$whitefate)
covid.white$whitefate_1 <- covid.white$whitefate/3
table(covid.white$whitefate_1)

table(covid.white$resentment)
covid.white$resentment_1 <- covid.white$resentment/16
table(covid.white$resentment_1)


pred_write1 <- lm(supportrelease_Rec3 ~  age_1 + woman + edu_1 + pid7_1, data = covid.white)
summary(pred_write1)

pred_write2 <- lm(supportrelease_Rec3 ~ age_1 + woman + edu_1 + pid7_1 + covidworry_1 + poordeserve_1 + whitefate_1 + resentment_1 , data = covid.white)
summary(pred_write2)


texreg(list(pred_write1, pred_write2), 
       custom.coef.names = c( 
         "(Intercept)",
         "Age",
         "Woman",
         "Education",
         "Partisanship",
         "Covid worry",
         "Poor deserve more",
         "White linked fate",
         "Racial resentment"
         
       ),
       digits = 3, custom.model.names = c("Model 1", "Model 2"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:predcontent", stars = c(0.05), caption = "Observational Analysis of Predictors of Writing a Pro-Release Message", float.pos = "t", single.row = TRUE)







